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IMAGING OF A SCATTERING MEDIUM USING 
THE EQUATION OF RADIATIVE TRANSFER 

This invention was made with U.S. Government support under contract number 
R01 AR46255-01, awarded by the National Institute of Health. The U.S. Government 
5 has certain rights in the invention. 

This application claims the benefit under 35 U.S.C. §120 of prior U.S. Provisional 
Patent Application Serial Nos. 60/177,779 filed January 24, 2000, entitled IMAGE 
RECONSTRUCTION IN OPTICAL TOMOGRAPHY. 

BACKGROUND OF THE INVENTION 

1 0 Field of the Invention 

The invention relates to the field of imaging in a scattering medium, and in 
particular to improved methods of reconstructing images of the spatial distribution of 
properties inside the scattering medium. 

Background Information 

15 Imaging of a scattering medium relates generally to an imaging modality for 

generating an image of the spatial distribution of properties (such as the absorption or 
scattering coefficients) inside a scattering medium through the detection of energy 
emerging from the medium. 

A typical system for imaging based on scattered energy measures, includes at 

20 least one source and detector on the surface of the medium for illuminating the medium 
and detecting emerging energy, respectively. The source illuminates the target medium 
with energy that is highly scattered in the medium so that the energy in general does not 
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travel along a straight line path through the medium, but rather, is scattered many times 
as it propagates through the medium. The detectors on the surface of the medium 
measure this scattered energy as it exits the target medium. Based on measurements of 
the scattered energy emerging from the target medium, a reconstructed image 
5 representation of the internal properties of the medium can be generated. These systems 
and methods are in contrast to projection imaging systems, such as x-ray , that measure 
and image the attenuation or absorption of energy traveling a straight line path between a 
the x-ray source and a detector. 

As will be appreciated, there are many instances where these techniques are 

10 highly desirable. For example, one flourishing application is in the field of optical 
tomography. In recent years optical methods using near-infrared energy (i.e., 
electromagnetic radiation with wavelengths in the range of -700— 1200 nanometers) 
have become increasingly important for noninvasive diagnostics in medicine. The ability 
to use near infrared (NIR) energy to probe the human body is of particular interest, 

15 because the propagation of NIR energy in tissue is exceptionally responsive to blood 
oxygenation levels and blood volumes, so that blood acts as a contrast agent. These 
attributes permit imaging of the vasculature, and thus provide great potential for detecting 
cardiovascular disease, tumors, brain hemorrhages, breast cancer, rheumatoid arthritis in 
finger joints, and the like. Furthermore, changes in the scattering properties may also be 

20 used as contrast agents. For example, rheumatoid arthritis affects the scattering 
coefficient of the synovial fluid, which is located inside the joints. 

Systems for making measurements on human subjects using sources and detectors 
as described above are well known and readily available. However, a major challenge 
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remains in the development of algorithms that transform these measurements into useful 
images on the target medium's internal properties. Since near-infrared energy is strongly 
scattered in tissue, standard backprojection techniques applied in X-ray tomography have 
been of limited success. Various other mathematical methods have been developed for 
5 solving the reconstruction problem in optical tomography. Some commonly applied 
schemes are modified backprojection methods, diffraction tomography, perturbation 
approaches, full matrix inversion techniques and elliptic system methods. Most of these 
methods employ model-based iterative imaging reconstruction (MOBIIR) techniques. 
These techniques employ a forward model that provides predictions of the detector 

10 readings based on an initial guess of the spatial distribution of properties inside the 
medium. The predicted detector readings are then compared with experimentally 
measured readings using an appropriately defined objective function <f>. The true 
distribution of internal properties of the target medium is then determined by iteratively 
updating the initial guess and performing new forward calculations with these updated 

15 internal properties until the predicted data agrees with the measured data from the 
detector readings. The final distribution of internal properties is then displayed or printed 
as an image. 

These known techniques, however, have several disadvantages. First, all of the 
known reconstruction schemes are generally based on the diffusion approximation to the 
20 equation of radiative transfer. The equation of radiative transport describes the 
propagation of photons through a scattering medium based in part on the internal 
properties of the target. The diffusion approximation to this equation, makes several 
assumptions that reduce the complexity of the equation making it easier and faster to 
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solve. However, the diffusion approximation is only valid for example: (1) for highly 
scattering media where the properties of the medium [i\ (the scattering coefficient) is 
much larger than (i a (the absorption coefficient), (2) for media that do not contain strong 
discontinuities in optical properties, such as very low scattering and absorbing regions 
5 ("void-like regions") embedded in highly scattering materials, and (3) for large source- 
detector pair separation, i.e., the separation between a source location and a detector 
location on the target medium. 

If the above conditions are not met, the diffusion approximation will fail to 
accurately describe energy propagation through the medium. Examples of highly 

10 absorbing regions in the body, where the scattering coefficient is not much larger than the 
absorption coefficient, are hematoma or liver tissue. Void like regions, those having low 
scattering and low absorbing areas, are present in areas such as the cerebrospinal fluid of 
the brain, the synovial fluid of human finger joints, or the amniotic fluid in the female 
uterus. For these cases the diffusion approximation fails and it is highly desirable to have 

15 a reconstruction code that is based on the theory of radiative transfer. 

Another disadvantage of the known techniques is that most of the currently 
available reconstruction schemes require the selection of a reference medium having 
properties nearly identical to the known properties of the target medium, so that there is 
only a small perturbation between the reference and target medium. Recently, increasing 

20 attention has been paid to gradient based iterative image reconstruction (GIIR) methods, 
which are a subclass of MOBIIR schemes. GIIR schemes overcome the limitations of 
small perturbations, by using information about the gradient of the objective function (the 
relative difference between predicted data for the reference and measured data from the 
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target) with respect to the distribution of optical properties in a minimization scheme. 
However, current techniques for calculating the gradient are impractical due to the 
complexity of the calculation and the corresponding time required to execute the 
calculation. 

5 For the foregoing reasons, there is a need for a reconstruction scheme that can 

efficiently and accurately reconstruct an image of the internal properties of a 
heterogeneous scattering medium, including but not limited to media that contain void- 
like regions and media that contain regions in which the absorption coefficient is not 
much smaller that the scattering coefficient. 

1 0 SUMMARY OF THE INVENTION 

The present method and system satisfies these needs by providing a gradient- 
based iterative reconstruction algorithm for imaging of scattering media using (1) the 
equation of radiative transfer as a forward model, and (2) an adjoint differentiation 
method for the fast and efficient gradient calculation used in the modification of the 

15 initial guess in the updating scheme. 

One embodiment of the present method and system provides a method for making 
an initial guess of the internal properties of a target medium, predicting the propagation 
of energy through the medium based on the initial guess, measuring the actual 
propagation of energy through the medium, updating the initial guess based on the 

20 predicted data and measured data. The predicted propagation of energy is calculated 

using the equation of radiative transfer and an initial guess as to what the spatial 

properties of the target might be. The measured propagation is obtained by directing 

energy into the target and measuring the energy emerging from the target. An objective 
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function is then generated that relates the predicted values to the measured values. The 
objective function gives an indication of how far the initial guess of the spatial properties 
of the target was from the actual spatial properties of the target. A gradient of the 
objective function is then generated using a method of adjoint differentiation. The 
5 gradient indicates how that initial guess should be adjusted to make the predicted data 
more closely match the measured data. This process is then repeated until the predicted 
and measured data are within an acceptable error. At this point the adjusted initial guess 
is representative of the spatial properties of the actual target and an image is generated 
therefrom. 

10 The foregoing specific objects and advantages of the invention are illustrative of 

those which can be achieved by the present invention and are not intended to be 
exhaustive or limiting of the possible advantages that can be realized. Thus, the objects 
and advantages of this invention will be apparent from the description herein or can be 
learned from practicing the invention, both as embodied herein or as modified in view of 

15 any variations which may be apparent to those skilled in the art. Accordingly, the present 
invention resides in the novel parts, constructions, arrangements, combinations and 
improvements herein shown and described. 

BRIEF DESCRIPTION OF THE DRAWING 

For a better understanding of the invention, together with the various features and 

20 advantages thereof, reference should be made to the following detailed description of the 

preferred embodiments and to the accompanying drawings wherein: 

FIG. 1 is a flow chart of the reconstruction process; 

FIG. 2 is a schematic of an exemplary measurement and imaging system; 
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FIG. 3 A is a computational graph of the forward model calculating the objective 
function; and 

FIG. 3B is a computational graph of the adjoint differentiation technique. 
DETAILED DESCRIPTION OF THE INVENTION 

5 I. Introduction 

As discussed above, imaging of a scattering medium refers generally to the 
methods and techniques of reconstructing an image of the internal properties of a 
scattering medium based on the transmission of energy into the medium and the 
measurement of the scattered energy emerging from the medium. 

10 The present method and system employs an iterative image reconstruction scheme 

that has three major elements. The following sections describe in detail how the three 
major elements of the present method and system work together to yield a reconstructed 
image of the internal properties of the scattering target medium. First, the solution of the 
forward model of the equation of radiative transfer is solved using an upwind difference 

15 scheme, even-parity finite-element formulation or the like; by way of example, the 
present invention is illustrated using the upwind scheme. This is followed with a detailed 
description of the second and third component, namely the objective function and the 
updating scheme, in which adjoint differentiation is used for the gradient calculation. 
Referring to FIG. 1, these elements are illustrated as the forward model 105, the objective 

20 function/analysis scheme 120 and the updating scheme 125. 

The forward model starts at step 100 with an initial guess [Xq = [|a s .o(r), M<a.o(r)] 
(scattering and absorption coefficients respectively) of the internal properties, known 
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energy source S(r s ) at the positions r s and given boundary conditions. Based on the initial 
guess \iq 9 energy source S(r s ) and the given boundary conditions the forward model (1) 
gives a numerical solution of the energy distribution T(r) inside the target scattering 
medium, and (2) predicts the energy radiance ^ on the boundary Q of the medium based 
5 on the equation of radiative transfer. Referring to step 110 of FIG. 1, the vector P, 
generated by the forward model, contains the predicted radiance value or derived 
parameter related thereto, at each detector, for each source detector pair, as a function of 
the properties |u. Measured data for the actual target to be imaged is then collected in 
step 115. The measured radiance at each detector, for each source detector pair is stored 
10 in the vector M. A typical imaging system for collecting the measured data is presented 
in detail below. 

In the analysis scheme, step 120, the predicted radiances ^((lo) are compared 
with the set of measured radiances M on the boundary Q of an actual target medium to 
define an objective function O. In steps 125 through 135, an optimization technique is 

15 used to minimize the objective function <I>, e.g., the normalized difference between the 
predicted and measured values. For this optimization the internal properties \x 0 are 
updated using the derivative dO/d// of the objective function with respect to the internal 
properties. In step 130, the gradient is computed by an adjoint differentiation technique, 
that takes advantage of a well-chosen numerical discretization of the forward model. 

20 Once the gradient is computed, the initial guess is modified in step 140 or 145, 

depending on whether an inner or outer iteration is being performed, and a new forward 
calculation is performed based on the new set of internal properties jlxo + A//. The 
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iteration process continues until the minimum of the objective function is reached within 
a specified error. At this point the predicted detector readings are identical to the 
measured detector readings within a given tolerance. The internal properties, //, of the 
target medium are then mapped into an image. 
5 For illustration purposes, the present system and method is described in further 

detail below with respect to an optical tomography system. However, it will be 
appreciated by those skilled in the art that the methodology of the present invention is 
applicable for any energy source (e.g., electromagnetic, acoustic, etc.), any scattering 
medium (e.g., body tissues, oceans, foggy atmospheres, geological strata, and various 

10 materials, etc.), and any source condition (e.g., time-independent, time-harmonic, time- 
resolved). Its applicability is dependent only on the presence of the phenomenology 
described herein, (i.e., radiative transfer being the principal mechanism of energy 
transport), which is expected in all cases where scattering occurs. Accordingly, this 
methodology can be extended to allow for new imaging approaches in a broad range of 

15 applications, including nondestructive testing, geophysical imaging, medical imaging, 
and surveillance technologies. 

II. Exemplary Measuring System 

There are many known imaging systems for collecting the measured data used in 
image reconstruction. A schematic illustration of an exemplary optical tomography 
20 system is shown in FIG. 2. This system includes a computer 202, energy source 204, a 
fiber switcher 208, an imaging head 210, detectors 212, a data acquisition board 214, 
source fibers 220 and detector fibers 222. 
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The energy source 204 provides optical energy, directed through the fiber 
switcher 208 to each of the plurality of source fibers 120 one at a time in series. The 
source fibers 220 are arranged around an imaging head 210 so that the energy is directed 
into the target medium 216 at a plurality of source locations around the target. 
5 The energy leaves the source fiber 220 at the imaging head 210 and enters the 

target medium 216 centered in the imaging head 210. The energy is scattered as it 
propagates through the target medium, emerging from the target medium at a plurality of 
locations. The emerging energy is collected by the detector fibers 222 arranged around 
the imaging head 210. The detected energy then travels through the detector fibers 222 to 

10 detectors 212 having energy measuring devices that generate a signal corresponding to 
the measurement. The data acquisition board 214 receives the measurement signal for 
delivery to the computer 202. 

This process is repeated for each source position so that a vector of measures are 
obtained for all of the detectors and source locations. This vector of measured data is the 

15 vector M referred to above. The vector of data M is stored by the computer 202 for use 
in image reconstruction and other analysis discussed below. 

III. Forward Model 

As discussed in the introduction above, the present method and system compares 
measured data M from the actual target with predicted data P for the target. The 
20 predicted data is obtained through the "forward" model. 

In the present method and system, the equation of radiative transfer is used for the 

forward model. Referring to FIG. 1, steps 100 through 110, the radiative transfer 

equation is used to predict the detector readings of energy emerging at one or more 
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detector locations on the scattering medium based on an initial guess of the properties of 
the medium. The equation of radiative transfer is derived by considering energy 
conservation in a small volume. The equation for time-dependent case can be written as: 

*^^ = S(r^ (la) 
5 and; for the time-independent case can be written as: 

<oW(r,<o) + (^ + \i s y¥(r 9 <Q)=S(r,G>) + \i s ^ n p^^Wr^d^ 9 , (lb) 
where r is the spatial position vector, co is a unit vector pointing in the direction of 
photon propagation, YfooM) and T(r 3 co) are the energy radiance in units of W cm" 2 sr" 1 . 
S(r ? G>,t) and S(r,co) are the source terms representing energy directed into a solid angle 
10 centered on oo in a unit volume at r, /u s is the scattering coefficient given in units of cm" 1 , 
ju a is the absorption coefficient given in units of cm" 1 and /?(g),g>') is the phase function 
that describes the probability that a photon with direction co 1 will be scattered into the 
direction co during a scattering event. 

The scattering phase function /?(©,gd ! ) can be described, for example, using the 
1 5 well known Henyey-Greenstein scattering function, 

/>(cos 0) = . (2) 

2(l + g 2 -2gcos0) 3/2 

where &is the angle between the two directions co and a) '. The parameter g is the 
anisotropy factor which is commonly used to characterizes the angular distribution of 
scattering. It can range from g=-l (total backward scattering), over g=0 (isotropic 
20 scattering), to g=l (strictly forward scattering). Biological tissues typically have a g- 
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factor between g=0.7 and g=0.98. Other scattering phase functions are also possible, for 
example derived from Mie theory. 

The integral of the radiance over all angles co at one point r yields the fluence rate 
O(r) for the time-independent case: 

5 0(r)= J Y(r,©) Jco (3) 

2% 

While it will be appreciated to those skilled in the art that a similar three- 
dimensional and/or time dependent equation can be generated, by way of example the 
remainder of the specification will focus on the two-dimensional, time-independent 
problem. 

10 Various computational techniques exist, such as those disclosed by Lewis E E, 

Miller W F, Computational Methods of neutron transport, New York, John Wiley & 
Sons, 1984, that numerically solve the equation of radiative transfer. Techniques 
commonly applied include the singular eigenfimction method, the method of spherical 
harmonics, the method of characteristics, the finite-element method, and the finite- 

15 difference discrete-ordinate method. A concise review of these techniques has been 
presented by Sanchez R, McCormick N J, A Review of neutron transport approximations. 
Nucl. Sci. Eng. 1982,80:481-535; and McCormick NJ, Inverse radiative transfer 
problems: a review, Nucl. Sci. Eng. 1992,112:185-198. The discrete-ordinates method is 
widely used with several different finite-difference approximations, such as the diamond- 

20 difference scheme, the weighted diamond-difference scheme, or the centered-difference 
scheme. 

Another means of solving the equation of radiative transfer is the upwind- 
difference scheme used in connection with the discrete-ordinates method to the equation 
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of radiative transfer disclosed by Sewell G., The numerical solution of ordinary and 
partial differential equations, San Diego, Academic Press, 1988. This method is a 
computationally efficient method for the calculation of the radiance and has the 
advantage that it easily supplies the required mathematical structure for the adjoint 
5 differentiation calculation. The adjoint differentiation method is crucial for obtaining the 
gradient of the objective function in an computationally efficient way, and thus obtaining 
an update of the initial guess of the properties of the target medium. The adjoint 
differentiation method and the objective function will be described in detail below in 
Section IV. 

10 Returning to the forward problem, to solve the equation of radiative transfer with 

an upwind-difference discrete-ordinates method, the angular and spatial variables have to 
be discretized. First, the integral term in equation lb is replaced by a quadrature formula 
that uses a finite set of K angular directions (Du with k= 1, .... K. This yields a set of K 
coupled ordinary differential equations for the angular-dependent radiance 

15 ^(r) = v F(r 3 G)k) in the directions 6%. The coupling term is the internal source term 

K / \ 

V> s X ctfc pip £.,(Q £/F(r 5 G)') . The parameter is a weighting factor that depends on the 
*'=i 

chosen quadrature formula. In this work the extended trapezoidal rule is employed. 

k 

®k V Vk( r ) + (^a +V s )V k (r) = S k (r) + ^ Z k ^ (4) 

Additionally, the spatial variable r needs to be discretized. The domain Q. is 
20 defined by a rectangular spatial mesh with M grid points on the x-coordinate and N grid 
points on the ^-coordinate. The distance between two adjacent grid points along the x- 
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axis is Ax and along the y-axis is Ay. The angular radiance at the grid point (ij) with 
position r = (x h yj) for a particular direction cok is represented by W^j = (xi, yj). The 
angular direction a>k can be expressed in cartesian coordinates with ^ = e x .0%= cos(^) 
and ?]k = ZytDk= sin(^t). The transport equation can now be written as: 

^ g^^O^+^A^^ Kp kk ^^(x h yj) (5) 

Finally, the spatial derivatives have to be replaced with a finite-difference 
scheme. The upwind-difference formula depends on the direction a>k of the angular- 
dependent radiance *F*. Thus, the set of all angular directions (Ok are subdivided into four 
quadrants to generate four different difference formulas for the radiance T^y, depending 
10 on the sign of ^ and % : 



I) >0,T, t > 0: f-S^,,,, = A "8,^ -^T^ 1 < 6a > 

5x Ax: dy Ay 

ii) 5* <o,m>0: 2_ „, 5 ^ ti j = , £_ H = m 

m) ^ > 0,n* < 0: ^ « 8^, „ = ^ ^ , |; « 8^ = (6c) 

15 IVUyfc<0,Tu<0 : — *8 = ^,-4^ 

> ^ dx x k -'- J Ax dy y K,J 

The time-independent radiative transfer equation, with the external and internal 
source terms on the right-hand-side, can now be written as: 





Ay 




J+l T l.7 




Ay 








Ay 
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Recasting the left-hand side of the preceding equation as a single operator acting 
upon ^k.u, produces; 

A + tu • 6,, + (ji fl + [i s )}¥ klJ = S KlJ + ^JX/V*^ , (8) 



from which it is apparent that the system of equations 8 corresponding to all K 
5 directions can be written as a single matrix equation: 

A*F= b; (9) 
proportional to one row of the matrix A and with an ordering of the radiance vector for a 
fixed k with, for example % k >0,iu>0, Y = 0P* )U ,T M , 2 , ... ^ k ,2^ k , 2 , 2 , ... 
■.■?¥kMJ>d produces: 

The resulting system of equations for the radiance W k jj for all grid points is 
solved for each ordinate index k by a Gauss-Seidel method. Accordingly, the matrix A is 
split into a diagonal part D, an upper-triangular part U, and a lower-triangular part L ? 
with A = L + D + U. The original matrix equation 9 is now be written as either: 
15 (L + D + U)Y = b (11a) 

or, (D + L)*F - -U¥ + b (lib) 

In this case an upper matrix U does not exist. Thus, for example for the case 
&>0, %>(): 

{* + % + ^ + ^ KlJ ~ * ^ " Ay ^ = Sk " + ^W***" (12) 

20 diagonal matrix D lower matrix L b 
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Now the iterative form with the iteration matrix (D+L) is expressed as 

(D + L)T Z + 1 = -UY Z + b , (13) 
and for the case ^ > 0, % > 0: 

5 {f^ + ta -f^Vu -^ z+ Vi = ^ (J +^|X/^V,,,(14) 

k,, +^,t w^*/ + :b vrVl - + 7t Xi ' z+lfc -- 1 l 

^ = J k ^—n 1 • (15) 

All and ^^y-i are a l rea( ly calculated at the current iteration step 

because of the vector ordering. For all other ordinates 0% besides ^ > 0, % > 0, the 
radiance vector Y has to be re-ordered to get the same matrix structure with 
10 (D + hy¥ = -IFF + b . The iteration process is repeated until the error norm 

E = I ^l^j ~~^ki\ at ^ e S^d P°i nt (y) i s smaller than a defined s. 

A significant improvement in convergence speed can be achieved by a slight 
modification to the Gauss-Seidel method. The successive over-relaxation (SOR) method 
uses a relaxation parameter p with \ <p < 2 in order to correct the solution Wj* 1 of the 

15 Gauss-Seidel iteration, now denoted as*?^. The updated value of the SOR is a 

linear combination of the Gauss-Seidel value 4^ and the previously computed value 

^ z Klj of the SOR using: 
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The boundary conditions are treated between each successive iteration steps. 
Because of the refraction index mismatch at the air-tissue interface, the outgoing radiance 
is partly reflected on the boundary and only a fraction of that energy escapes the medium. 
The internally reflected energy contributes further to the photon propagation inside the 
5 medium, while the transmitted energy enters the detectors. Using Fresners formula, the 
transmissivity T and reflectivity R are calculated at the boundary grid points for each 
ordinate CDk and for the given refractive indices. The reflectivity R and the transmissivity 
Tare defined as 

R _ Sln2 fatnm, ~ <Q+ Sln2 trans + <Q ^ 

10 T = l-R (18) 

The angles © /raBJ , which pertain to radiances escaping the object, are determined 
by Snell's law (n object sin© IB = w^sin©^) given the angle © w of the radiance, which hits 
the boundary inside the object. The angle to re/ = -® m is just the angle of the reflected 

radiance on the boundary inside the object. The transmitted and the reflected radiance 
15 are calculated with: 

A trans 1 1 in ' V 1 y ) 

Xe f = R-X- (20) 

The fluence on the boundary at the grid point (ij) 9 which enters the detector, is 

calculated for a given detector aperture AP using the transmitted radiance, 

k=k 2 

20 *W = 1^W>^ * 2 a k W kiU =O hJ . (21) 

AP k=k x 
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The weighting factors are given by the extended trapezoidal rule, [84] which 
provides the quadrature formula in this work. 

These fluence values o & on the boundary provide the vector of predicted values 

shown in step 110 of FIG. 1 used for the reconstruction algorithm. 

5 IV. Objective Function/Analysis Scheme and Updating Scheme 

The first components of the present method and system, i.e., the forward model 
used to generate the predicted energy emerging from the target, and the collection of 
measured energy emerging from the target, have been described in detail above. The 
following sections describe the second and third components of the method and system, 

10 illustrated in steps 120 through 135 of FIG. 1, namely, the objective function and the 
updating scheme, in which adjoint differentiation is used for the gradient calculation. 
These second and third components are used to modify the initial guess of the internal 
properties of the medium in an iterative process. 

1. Objective Function/Analysis Scheme 

15 Optical tomography can be viewed as an optimization problem with a non-linear 

objective function. Referring to FIG. 1, step 120, the discrepancy between the actual 
measurements, represented by the vector M, and the predictions, given by the vector P, 
for m source-detector pairs is defined by a scalar-valued objective function <t>(T| The 
predictions P are mapped onto a scalar <p using the % 2 - error norm: 

20 0 : 9T -» 9? (22a) 

i m « 

P^9=TS(^-M / y (22b) 
18 
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The predictions P are determined for all m source-detector positions by evaluating 
the forward model F 9 given the spatial distribution of internal properties //(see Part I): 

F:9T-»9T (23a) 
\i h> P([i) (23b) 
5 The vector // is of length n =2 IJ, and contains the internal properties ju s and 

The parameters / and / denote the number of grid points of a finite-difference grid along 
the x-axis and y -axis, respectively. 

The goal in image reconstruction is to determine a vector // that minimizes the 
objective function d>. As already mentioned in the introduction, gradient-based 
10 optimization techniques, such as the conjugate-gradient method, employ the gradient 
V //]>(//) to calculate a sequence of points /4>, //u ...,/// that lead to ever-improving 
values of <E>. This iterative process is terminated when |<E>(/ft+7) - €>(///) | becomes smaller 
than a predefined value s. A crucial task within this process is the computationally 
efficient calculation of the gradient V^O^). The following section describes in detail 
15 how this can be done using an adjoint differentiation scheme. 

2. Gradient Calculation With Adjoint Differentiation 
In adjoint differentiation schemes the numerical code that calculates the objective 
function O is directly differentiated to obtain the gradient V A <D(/4), FIG. 1, step 130. To 
apply this method one has to decompose the objective function into a series of elementary 
20 differentiable function steps. By systematically applying the chain rule of differentiation 
to each of these elementary steps in the reverse direction of the forward model 
computation, a value for the gradient is obtained. 
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The parameter G> = 0(P(//)) can be decomposed into sub-functions F* in the 
following way: 

«(^))=4.(^(f-(f- ! (..(fXF , ^)»^» p4) 

: = (O o F z (n) o F z_1 (n) o F z ~ 2 (ji) c ... o F 2 (p) o F 1 )^). 

The sub-functions i 5 * are defined by the SOR-method that is used to solve the 
5 forward model (see Part I). The SOR-method is an iterative approach and the z-th 
iteration yields the intermediate result The radiance vector ^is of length p = 

KxIxJ with£ e [l,^], / € [l,/], and j The detector readings P(//) are the angular- 

dependent radiances at the last iteration step Z at detector positions (ij) on the 

boundary. The computational graph of the forward model is depicted in FIG 3 A. Starting 
10 with the internal properties ju as the input variables, the sub-function F l produces the 
intermediate result and output variable ¥*. The output variables ¥ of F* and the internal 
properties ju always serve as input variables for the next sub-function for all 
subsequent steps of the decomposition. The step F z calculates the predictions P, which 
become the input to the final step of the objective function <D determination, which is the 
15 calculation of the scalar q>. 

The sub-functions F map the intermediate variables and constant input 

values of ju onto the intermediate result T z = F z {^f z " \ \ij for all iteration steps of the 
transport forward model: 
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F'-M'xW -»9t' (25a) 
h-^ 2 , (25b) 



For the ordinates with ^ > 0, r\j c >0 this mapping is given explicitly by: 



= 0-<D>F* Vr 1 £=1 — TT 1 ~ (26) 

For the first iteration step z=l, only maps the input variables ju : 



F* : 9T -> 3l p (27a) 
H i-> (27b) 
and produces: 



^ V„ = (o . (28) 

Equations 26 and 28 are the smallest computational units in the forward code and 
play an important role in the adjoint differentiation step, which is discussed below. The 
gradient V^O of the objective function is obtained by differentiating equation 24 with 

respect to the internal properties //. The total derivative V^O can be obtained by 

systematically applying the chain rule of differentiation to equation 22. The total 
derivative V^O is a composition of derivatives of all intermediate steps of the forward 
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model due to the explicit dependence of the intermediate results on the internal 
properties. Starting from the last step of the forward code, in which the objective function 
O is calculated given the predicted detector readings P = x ¥ z , the total derivative is given 
by the chain rule as: 



JO 



f do y dv z 



d\i) Kd¥ z ) dp. + UuJ 



' 5® 



(29) 



The derivative 



d\i 



is obtained by again applying the chain rule of 



differentiation: 



dV z &¥ z di> z ~ x &¥ z 
d\i ~ d¥ z ~ l d\i + dp. 



(30) 



Equations 29 and 30 yield together: 



10 



'doY _(&!>Y fflT z dSV z ~ x ( c®Yd¥ z 



(31) 



In the next step the total derivative 



■2-1 



is replaced again. This process is 



repeated for all total derivatives down to the first step — — = — - , to obtain: 

l^J "1W Z J dy ZA '"dv l diirKff¥ z J dv 2 - 1 "^^ 2 diir m " + V<ff¥ z ) an J + UnJ 



15 



or, using the short hand notation 

[r doY fflF z ^1 9(fl>oY z o...o^) r (dd>X 



(33) 



rewriting equation 33 produces: 
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(d^YdV* ( d&Y 8V Z (d®" T 



dp) W 1 ^ dp + W 2 J dp + "' + {cw) dp + "' + la¥ z J dp + {dp 

(35) 



(34) 





T 


u 










\dp) 








+ 





or, V^O:= 



which contains three distinct terms. The last term | — 

. d\ij 



is zero, because O is not an 



explicit function of the internal properties. The middle term — — can be calculated from 

op- 
equation 26 of the forward model. The derivatives with respect to ju a and ju s are: 



K 



[Ax V ^ 



CO" 



- -co 



(36a) 



(36b) 



At this point the adjoint differentiation technique has not yet been used, since the process 
10 has not stepped backward through the forward code. This procedure comes into play in 

fdd>Y 

the calculation of the first term J in equation 35. This can be best understood 

while looking at the computational graph of the adjoint differentiation technique in FIG 
3B. Starting with the last step (output) of the forward code, which is the calculation of 
the objective function given the predictions P = i^, O is differentiated with respect to 
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7 ( &S?Y 

*F and multiply it with ^— j = 1 . The result is difference between the measured and 



predicted data, 



_5£ 



= (t z -mJ. 



(37) 



This is the input to our adjoint model, which will eventually provide the gradient 

( SO V 

5 V M <D. More specifically, \^ is calculated by continuing to step backward through 



the forward code and is given by: 

f do Y 



(dv z Y(d®Y 

{d*¥ z - l J ^ Z J • 



(38) 



(doY 

The remaining derivatives j °f a ^ intermediate steps in equation 35 are 

( d<D Y 

computed using the previously calculated derivatives j , such that 



(d&) T _ a(<Do,..o^ z+1 X x P z y 



10 



dv 2 



aT 2 



(39) 



This step, in which (~~~) is calculated from > constitutes the adjoint 



differentiation step in our code. The matrix 



is calculated by analytically 
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differentiating the (z+lj-th SOR-iteration step, given in equation 26, with respect to *P. 
We get, for example in the case of ordinates with \ k > 0,r\ k > 0 : 

gy... , . = (i + © J n — : — ; n ( 4 °) 



with 8.^=8^ and 8, ={J j£ 



5 Here, the approximations k^il±. s lt S and Jk ^ +1 ^ s ik 8 

can be made for the relevant terms on the right hand side of equation 40. 

As can be seen, the gradient of the objective function is calculated by stepping 
backward through all previously calculated iteration steps of the forward model without 
solving an entirely new numerical adjoint equation of radiative transfer. Furthermore, the 
10 particular underlying physical system does not have to be known, because the derivative 
is computed directly from the elementary steps of the forward model code (equations 26 
and 28). 

3. Gradient-Based Optimization 

Once the gradient V^O^) for a point jUq is obtained, a one-dimensional line 

15 minimization along the given gradient direction is performed until a minimum, <E>(/4ni n ) is 
found, FIG. 1, step 135 (inner iteration). Various techniques can be applied to perform 
such one-dimensional minimizations. Exemplary techniques are disclosed in Nash S G, 
Sofer A. Linear and nonlinear programming, McGraw-Hill, New York, 1996, and Press 
W H, Teukolsky S A, Vetterling W T, and Flannery B P. Numerical Recipes in C. 
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Cambridge University Press, New York, (1992). One approach is to start with three 
points, the initial guess jUo, ju\ = juo + Aju , and ju% = jtio +2Aju chosen along the direction 
of the gradient. After calculating the objective function at these three points, the largest 
calculated value is neglected, and a new guess /& is determined using the golden section 
5 rule until a minimum is bracketed. At that time a parabola is fit through these three 
points. The smallest value of the objective function on this parabola is determined to be 
the minimum. Once the point (guess) is found for which <£(//min,i) is minimal along 
the one-dimensional subspace, a new gradient calculation V / /X>(/4 n i tl( i) is performed at the 
minimum /4nm,i and a new direction is chosen in a conjugate-gradient framework. For 

10 this a Polak-Riberie conjugate-gradient scheme can be employed (outer iterations.) 

After completing several inner and outer iterations a final point /fei n ,fmai is found 
for which ^>{ju^ m ^ m2i \) is smaller than for all other points ju^ before. The final 
reconstruction image is than given by /w,fmai ? wherein //mi n ,fmai represents the spatial 
distribution of the optical properties in the target medium. The image can be represented 

1 5 by any known means, such as on a computer monitor or printed as a hard copy on paper, 
wherein the property values of the medium may be represented as shades of gray or 
colors. 

Although illustrative embodiments have been described herein in detail, those 
skilled in the art will appreciate that variations may be made without departing from the 
20 spirit and scope of this invention. Moreover, unless otherwise specifically stated, the 
terms and expressions used herein are terms of description and not terms of limitation, 
and are not intended to exclude any equivalents of the system and methods set forth in the 
following claims. 
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